// LAB DIRECTIVES

typedef long int lint;
typedef long unsigned int luint;


#define LOGTRICK
#define STUCKPARTICLE_FIX
#define BASISWIND_FIX
#define MASS

#define TRHO float

#define FLOAT

#define DG_DTMAX_CELL
#define DG_ADD_STD

#define DV_DH2 // DV_DH (default)

#define ICS_FIND_VEL "ICs::find_vel.h"
#define ICS_CALC_ORB "ICs::calc_orb_sse.h"
//#define ICS_CALC_ORB "ICs::calc_orb.h"
//#define ICS_FIND_VEL "ICs::find_vel_future_nosync_fgrid_Lzmax_7points.h"


#include "ICs.h"
		
class MyICs: public ICs {
	public:
		
		int argc; 
		char **argv;
		
		
		MyICs() : argc(0), argv(0) {
			printf("MyICs::MyICs()\n");
		}
		
		MyICs(int _argc, char *_argv[]) : 
		argc(_argc), argv(_argv)
		{
			printf("MyICs::MyICs(int _argc, char *_argv[])\n");
		}

		void run() {
			
			make_halo();
			//make_disk();
 			//make_gf1_ICs();
			
			//cals_and_save_relden("makeICs.B.halo-disk.1");
			//cals_and_save_zmean("makeICs.B.halo-disk.5");
		}
		
		
		void make_halo() {

			
			if (!ag.load("ag.halo.dump..dat")) {

				ag.nphi=2;
				ag.nsign =4;
				ag.G = G;
				ag.softening = 0.1;
				ag.ErrTolTheta = 0.35;
				
				ag.r1 = 1;
				ag.nr = 16*100;
				ag.r2 = pow(2,16);
				
				ag.z1 = 1;
				ag.nz = 16*100;				
				ag.z2 = pow(2,16);

				
				dg.r1 = 1;
				dg.nr = 16*8;
				dg.r2 = pow(2,16); 
	
				dg.z1 = 1;
				dg.nz = 16*8;
				dg.z2 = pow(2,16);
				
				ns =  256; // steps per orbit
				np =   64; // periods per orbit		

				halo_rmax = ag.r2;
				halo_n = 10000;
			
				disk_mu = 0.04;
				disk_r0 = 2.5;
				disk_rmax = 100;
				disk_z0 = .5;
				disk_zmax = 100;
				disk_n = 0;
			
				bulge_mu = 0.0;
				bulge_r0 = .1;
				bulge_rmax = 1000;		
				bulge_n = 0;
				disk_mu = 0.00;
				
				fLz = 0;
				dvk_mu = 0.0025;
				dvh_mu = 0.0025;

				
				read_parameters(argc, argv);
				
				update_derivatives();
				
				//print();
				
				/*
				halo.n = 1e7;
				halo.p = new P[halo.n];
				halo.print();
				halo.set_pos(halo.p, halo.n);
				halo.save_pos("halo.pos.dat");
				ag.calc(halo.p, halo.n);
				halo.free();
				/**/

				
				
				for (int iz=0; iz<ag.nz+1; iz++) 
					for (int ir=0; ir<ag.nr+1; ir++) 
					{

						
						INT ri = ag.r1i + ir*ag.dri;
						INT zi = ag.z1i + iz*ag.dzi;

						REAL r = -ag.r1 + *(REAL*)&ri;
						REAL z = -ag.z1 + *(REAL*)&zi;
						
						//printf("%i %i: %g %g\n", iz, ir, r, z);
						
						double R = sqrt(SQR(r) + SQR(z));
						double A = halo.acc(R);
						
						AG::A &ai = ag.a[iz*(ag.nr+1) + ir];
						double A_over_R = A/R;
						
						ai.r = r*A_over_R; 
						ai.z = z*A_over_R;
						ai.pot = halo.pot(R);
						
						if (ir==0) ai.r=0;
						if (iz==0) ai.z=0;
				}
				
				/*
				Disk disk = this->disk;
				disk.n = 1e7;
				disk.p = new P[disk.n];
				disk.print();
				
				disk.set_pos(disk.p, disk.n);
				ag.calc(disk.p, disk.n);
				disk.free();
				*/
				ag.save("ag.halo.dump.dat");
				
			}
			
			//ag.save_txt("ag.halo.txt");	
			//AG::A a = ag.get_acc(1.0, 0.01);
			//a.print();
			

			// rho0
			m0 = new REAL[dg.n];
			
			printf("%i\n", dg.n);
			
			dg.calc(halo, m0);
			//dg.save_txt("rho0.halo.a.txt", rho0);
			//dg.calc(halo.p, halo.n, rho0);
			//dg.save_txt("rho0.disk.p.txt", rho0);
			
			//save_arr("m0.halo.dat", m0, dg.n);

			
				
			//exit(0);
			
			
			//exit(0);
			
			/*
			double *rho0_2 = new double[dg.size()];
			
			disk.G = G;
			disk.n = 20000000;
			disk.p = new P[disk.n];
			//disk.rmax = 10;
			disk.set_pos(disk.p, disk.n);
			disk.print();
			//disk.set_vel_uniform(disk.p, disk.n);
			dg.calc(disk.p, disk.n, rho0_2);
			dg.save_txt("rho0.disk.p.txt", rho0_2);
			*/
			
			trackF = false;
			track_freq = 1;
			tracks_freq = 500;
			
			

			
			halo.p = new P[halo.n];
			halo.set_pos(halo.p, halo.n);
			halo.set_pot(halo.p, halo.n);
			halo.set_vel(halo.p, halo.n);			
			halo.save_vel("halo.vel.dat");

			
			

			
			p = halo.p;
			n = halo.n;
			
			wn = sg.n*vg.n;
			w  = new float[wn];
			wa  = new float[wn];
			
			

			
			calc_w(wa);
			save_w(wa, "w.analytic.a.dat");
			
			
			halo.set_vesc(halo.p, halo.n, ag);
			halo.set_vel_uniform(halo.p, halo.n, ag);
			
			//halo.load_vel("tmp.halo.only.dh_min0.1dvmu0.01DH2n50k/vel/vel.00013.dat");
			//halo.load_vel("tmp/vel/vel.00009.dat");
			//save_gadget1("snap.halo.gf1");
			
			
			//halo.print();
			
			
			

			
			find_vel();
		
		}		

		
		void resume() {
			
			FileListFinder flf;
				
			flf.find(curfolder, "dump");
			
			if (flf.list.size() == 0) {
				printf("resume() :: ERROR :: no dump file\n");
				exit(0);
			}
			
			load( flf.list[0].c_str() );
			resumeF = true;

			m0 = new REAL[dg.size()];
			dg.calc(halo, m0);
			save_arr("rho0.dat", m0, dg.n);
						
			w  = new float[wn];
			nsteps = 10000000;
			find_vel();

		}

		

	
		void free_arrays() {

			delete [] m0;
			//delete [] p;
		}
	
	
	void read_parameters(int argc, char *argv[]) {
	
		if (argc==0) return;
		
		class ATO { public: // ...just to simplify what is bellow
			char **argv; string arg;
			ATO(char **_argv): argv(_argv) {}
			double f(int& i) { arg+=string(argv[i])+string(argv[i+1])+string(".");  return atof(argv[++i]); }
			int i(int& i) { arg+=string(argv[i])+string(argv[i+1])+string("."); return atoi(argv[++i]); }
		} ato(argv);
		
		for (int i=1;i<argc;i++)
				  if (!strcmp(argv[i],"np"))  np = ato.f(i); 
			else if (!strcmp(argv[i],"ns"))	ns = ato.f(i); 
			else if (!strcmp(argv[i],"dg.nr"))	dg.nr = ato.i(i);
			else if (!strcmp(argv[i],"dg.nz"))	dg.nz = ato.i(i);
			else if (!strcmp(argv[i],"dg.n")) 	dg.nr = dg.nz = ato.i(i);
			else if (!strcmp(argv[i],"dg.rz1")) dg.r1 = dg.z1 = ato.f(i);
			else if (!strcmp(argv[i],"dg.rz2")) dg.r2 = dg.z2 = ato.f(i);
			else if (!strcmp(argv[i],"n"))		halo_n = ato.i(i);
			else if (!strcmp(argv[i],"h"))		h = ato.f(i); 
			else if (!strcmp(argv[i],"dvk_mu"))	dvk_mu = ato.f(i); 
			else if (!strcmp(argv[i],"dvh_mu")) dvh_mu = ato.f(i); 
			else if (!strcmp(argv[i],"stype")) 	stype = ato.i(i); 
			else if (!strcmp(argv[i],"nsteps")) nsteps = atoi(argv[++i]); 
			else if (!strcmp(argv[i],"resume")) resumeF = true;	
			else if (!strcmp(argv[i],"ssfolder")) 	sprintf(ssfolder, "%s", argv[++i]);
			else if (!strcmp(argv[i],"curfolder")) sprintf(curfolder, "%s/", argv[++i]);
			
			
			if (ato.arg.length() && string(curfolder)==string("./"))  
				sprintf(curfolder, "%s/", ato.arg.substr(0, ato.arg.length()-1).c_str());
			
	}
};



int main(int argc, char *argv[]) {
	
	//std::cout << 8/-3 << std::endl;
	//exit(0);
	
	MyICs ics(argc, argv);
	
	ics.read_parameters(argc, argv);
	ics.resumeF ? ics.resume() : ics.run();
	
	return 0;
	
}